subroutine create_object(lid,px,py,pz,ierr,dx,dy,dz)
  use param
  use event_param
  use mat_param,only:half_a0
  use table_param,only:sim_obj
  use table_param,only:n_sim_obj
  use table_param,only:event_table
  use table_param,only:n_event_table
  use table_param,only:alloc_event_n_in_sim_obj
  use table_param,only:event_profile
  use table_param,only:n_event_profile
  use table_param,only:alloc_event_n_in_profile
  use table_param,only:spont_event_table
  use table_param,only:n_spont_event_table
  use box_param,only:vbox
  use box_param,only:vbox_map
  use box_param,only:cell_width
  use box_param,only:ch_width
  use label_data_param,only:obj_type
  use react_data_param,only:react_type
  use react_data_param,only:react_map
  use profile_map_param,only:profile_map_diff
  use profile_map_param,only:profile_map_diss
  use profile_map_param,only:profile_map_rctn
  use profile_map_param,only:profile_map_sink
  use geom_data_param,only:geom_obj
  use geom_data_param,only:n_geom_obj
  use sink_data_param,only:sink_map
  use sink_data_param,only:sink_type
  use count_param,only:count_obj
  use count_param,only:total_count_obj
  implicit none

  interface
     function reallocInt(p,n)
       integer,pointer,dimension(:) :: reallocInt,p
       integer,intent(in) :: n
     end function reallocInt
  end interface

  integer,intent(in) :: lid,ierr
  real(dbl),intent(in) :: px,py,pz,dx,dy,dz

  logical :: removed
  integer :: i,j,k,n
  integer :: dir,cur,event_id
  integer :: sid,bid,eid,pid,nid,rid,gid
  real(dbl) :: hx,hy,hz,xi,dist2,npos
  integer :: nx,ny,nz

  ! Add to the sim_obj
  n_sim_obj=n_sim_obj+1
  call alloc_sim_obj(n_sim_obj,sid)

  ! index
  sim_obj(sid)%index=sid

  ! object type
  sim_obj(sid)%object_type=lid

  ! position
  sim_obj(sid)%pos(1)=px
  sim_obj(sid)%pos(2)=py
  sim_obj(sid)%pos(3)=pz

  ! alignment
  if (ierr.eq.0) then
     sim_obj(sid)%vec(1)=dx
     sim_obj(sid)%vec(2)=dy
     sim_obj(sid)%vec(3)=dz
     sim_obj(sid)%dr(1)=dx*half_a0
     sim_obj(sid)%dr(2)=dy*half_a0
     sim_obj(sid)%dr(3)=dz*half_a0
  else
     call random_number(xi)
     dir=ceiling(xi*8.0)
     select case(dir)
     case(1)
        sim_obj(sid)%vec(1)=1
        sim_obj(sid)%vec(2)=1
        sim_obj(sid)%vec(3)=1
        sim_obj(sid)%dr(1)=half_a0
        sim_obj(sid)%dr(2)=half_a0
        sim_obj(sid)%dr(3)=half_a0
     case(2)
        sim_obj(sid)%vec(1)=-1
        sim_obj(sid)%vec(2)=1
        sim_obj(sid)%vec(3)=1
        sim_obj(sid)%dr(1)=-half_a0
        sim_obj(sid)%dr(2)=half_a0
        sim_obj(sid)%dr(3)=half_a0
     case(3)
        sim_obj(sid)%vec(1)=-1
        sim_obj(sid)%vec(2)=-1
        sim_obj(sid)%vec(3)=1
        sim_obj(sid)%dr(1)=-half_a0
        sim_obj(sid)%dr(2)=-half_a0
        sim_obj(sid)%dr(3)=half_a0
     case(4)
        sim_obj(sid)%vec(1)=1
        sim_obj(sid)%vec(2)=-1
        sim_obj(sid)%vec(3)=1
        sim_obj(sid)%dr(1)=half_a0
        sim_obj(sid)%dr(2)=-half_a0
        sim_obj(sid)%dr(3)=half_a0
     case(5)
        sim_obj(sid)%vec(1)=1
        sim_obj(sid)%vec(2)=1
        sim_obj(sid)%vec(3)=-1
        sim_obj(sid)%dr(1)=half_a0
        sim_obj(sid)%dr(2)=half_a0
        sim_obj(sid)%dr(3)=-half_a0
     case(6)
        sim_obj(sid)%vec(1)=-1
        sim_obj(sid)%vec(2)=1
        sim_obj(sid)%vec(3)=-1
        sim_obj(sid)%dr(1)=-half_a0
        sim_obj(sid)%dr(2)=half_a0
        sim_obj(sid)%dr(3)=-half_a0
     case(7)
        sim_obj(sid)%vec(1)=-1
        sim_obj(sid)%vec(2)=-1
        sim_obj(sid)%vec(3)=-1
        sim_obj(sid)%dr(1)=-half_a0
        sim_obj(sid)%dr(2)=-half_a0
        sim_obj(sid)%dr(3)=-half_a0
     case(8)
        sim_obj(sid)%vec(1)=1
        sim_obj(sid)%vec(2)=-1
        sim_obj(sid)%vec(3)=-1
        sim_obj(sid)%dr(1)=half_a0
        sim_obj(sid)%dr(2)=-half_a0
        sim_obj(sid)%dr(3)=-half_a0
     end select
  end if

  ! index for vbox
  nx=ceiling(px/cell_width)
  ny=ceiling(py/cell_width)
  nz=ceiling(pz/cell_width)
  bid=vbox_map(nx,ny,nz)
  sim_obj(sid)%vbox_index=bid

  ! initialize linked list and .removed.
  sim_obj(sid)%next=0
  sim_obj(sid)%prev=0
  sim_obj(sid)%removed=.false.

  ! save the index of event_table_list in sim_obj
  event_id=sim_obj(sid)%n_event

  ! ---- Treat the associated events ---- !

  ! event type 1 : diffusion
  if (obj_type(lid)%does_move) then

     ! add to the event table
     n_event_table=n_event_table+1
     call alloc_event_table(n_event_table,eid)

     ! record the index of event table in the sim_obj
     event_id=event_id+1
     sim_obj(sid)%n_event=event_id
     if (event_id.gt.sim_obj(sid)%n_alloc_event) then
        sim_obj(sid)%n_alloc_event=&
             sim_obj(sid)%n_alloc_event+alloc_event_n_in_sim_obj
        sim_obj(sid)%event_table_list=>&
             reallocInt(sim_obj(sid)%event_table_list,&
             sim_obj(sid)%n_alloc_event)
     end if
     sim_obj(sid)%event_table_list(event_id)=eid

     ! Set the event table
     event_table(eid)%obj1=sid
     event_table(eid)%obj_seq1=event_id
     event_table(eid)%event_type=diff_event
     event_table(eid)%removed=.false.
     event_table(eid)%is_spontaneous=.false.

     ! Set the event profile
     if (obj_type(lid)%diff_const_rate.eq.0) then
        ! if the event is spontaneous,
        ! add to the spont_event_table
        event_table(eid)%is_spontaneous=.true.
        n_spont_event_table=n_spont_event_table+1
        spont_event_table(n_spont_event_table)=eid
        event_table(eid)%spont_table_index=&
             n_spont_event_table
     else
        ! if the event is probabilistic,
        ! add to the event_profile
        if (profile_map_diff(lid).eq.0) then
           ! Create a new profile
           n_event_profile=n_event_profile+1
           call alloc_profile_table(n_event_profile,pid)
           profile_map_diff(lid)=pid
           event_profile(pid)%event_type=diff_event
           event_profile(pid)%rate_const=&
                obj_type(lid)%diff_const_rate
           event_profile(pid)%index=lid
           event_profile(pid)%n_event=1
           if (event_profile(pid)%n_alloc_event.eq.0) then
              event_profile(pid)%n_alloc_event=&
                   alloc_event_n_in_profile
              event_profile(pid)%event_list=>&
                   reallocInt(event_profile(pid)%event_list,&
                   event_profile(pid)%n_alloc_event)
           end if
           event_profile(pid)%event_list(1)=eid
           event_table(eid)%profile_table_index=pid
           event_table(eid)%profile_table_seq=1
        else
           ! Add to the exisiting profile
           pid=profile_map_diff(lid)
           n=event_profile(pid)%n_event
           n=n+1
           event_profile(pid)%n_event=n
           if (n.gt.event_profile(pid)%n_alloc_event) then
              event_profile(pid)%n_alloc_event=&
                   event_profile(pid)%n_alloc_event+&
                   alloc_event_n_in_profile
              event_profile(pid)%event_list=>&
                   reallocInt(event_profile(pid)%event_list,&
                   event_profile(pid)%n_alloc_event)
           end if
           event_profile(pid)%event_list(n)=eid
           event_table(eid)%profile_table_index=pid
           event_table(eid)%profile_table_seq=n
        end if
     end if
  end if

  ! event type 2 : dissociation
  do i=1,obj_type(lid)%diss_type

     ! add to the event table
     n_event_table=n_event_table+1

     ! record the index of event table in the sim_obj
     call alloc_event_table(n_event_table,eid)
     event_id=event_id+1
     sim_obj(sid)%n_event=event_id
     if (event_id.gt.sim_obj(sid)%n_alloc_event) then
        sim_obj(sid)%n_alloc_event=&
             sim_obj(sid)%n_alloc_event+alloc_event_n_in_sim_obj
        sim_obj(sid)%event_table_list=>&
             reallocInt(sim_obj(sid)%event_table_list,&
             sim_obj(sid)%n_alloc_event)
     end if
     sim_obj(sid)%event_table_list(event_id)=eid

     ! Set the event table
     event_table(eid)%obj1=sid
     event_table(eid)%obj_seq1=event_id
     event_table(eid)%event_type=diss_event
     event_table(eid)%sqnc=i
     event_table(eid)%removed=.false.
     event_table(eid)%is_spontaneous=.false.

     ! Set the event profile
     if (obj_type(lid)%diss_const_rate(i).eq.0) then
        ! if the event is spontaneous,
        ! add to the spont_event_table
        event_table(eid)%is_spontaneous=.true.
        n_spont_event_table=n_spont_event_table+1
        spont_event_table(n_spont_event_table)=eid
        event_table(eid)%spont_table_index=&
             n_spont_event_table
     else
        ! if the event is probabilistic,
        ! add to the event_profile
        if (profile_map_diss(lid,i).eq.0) then
           ! Create a new profile
           n_event_profile=n_event_profile+1
           call alloc_profile_table(n_event_profile,pid)
           profile_map_diss(lid,i)=pid
           event_profile(pid)%event_type=diss_event
           event_profile(pid)%rate_const=&
                obj_type(lid)%diss_const_rate(i)
           event_profile(pid)%index=lid
           event_profile(pid)%sqnc=i
           event_profile(pid)%n_event=1
           if (event_profile(pid)%n_alloc_event.eq.0) then
              event_profile(pid)%n_alloc_event=&
                   alloc_event_n_in_profile
              event_profile(pid)%event_list=>&
                   reallocInt(event_profile(pid)%event_list,&
                   event_profile(pid)%n_alloc_event)
           end if
           event_profile(pid)%event_list(1)=eid
           event_table(eid)%profile_table_index=pid
           event_table(eid)%profile_table_seq=1
        else
           ! Add to the existing profile
           pid=profile_map_diss(lid,i)
           n=event_profile(pid)%n_event
           n=n+1
           event_profile(pid)%n_event=n
           if (n.gt.event_profile(pid)%n_alloc_event) then
              event_profile(pid)%n_alloc_event=&
                   event_profile(pid)%n_alloc_event+&
                   alloc_event_n_in_profile
              event_profile(pid)%event_list=>&
                   reallocInt(event_profile(pid)%event_list,&
                   event_profile(pid)%n_alloc_event)
           end if
           event_profile(pid)%event_list(n)=eid
           event_table(eid)%profile_table_index=pid
           event_table(eid)%profile_table_seq=n
        end if
     end if
  end do

  ! event type 3 : reaction
  do i=1,vbox(bid)%n_neighbor

     ! Sweeping over the neighbor cells
     nid=vbox(bid)%neighbor_list(i)
     cur=vbox(nid)%first

     ! Loop over the objects in the cell
     do j=1,vbox(nid)%n_obj

        ! Check if the reaction is defined in the data sheet.
        rid=react_map(lid,sim_obj(cur)%object_type)
        if (rid.eq.0) go to 100

        ! Calculate the distance between objects
        hx=sim_obj(cur)%pos(1)+vbox(bid)%dx(i)-px
        hy=sim_obj(cur)%pos(2)+vbox(bid)%dy(i)-py
        hz=sim_obj(cur)%pos(3)+vbox(bid)%dz(i)-pz
        dist2=hx*hx+hy*hy+hz*hz

        ! Loop over the possible reaction types
        do k=1,react_type(rid)%n_data

           ! Compare the distance and capture radius
           if (dist2.gt.react_type(rid)%react_rd2(k)) cycle

           ! Determine the index of the event table
           n_event_table=n_event_table+1
           call alloc_event_table(n_event_table,eid)

           ! set the object 1
           event_id=event_id+1
           sim_obj(sid)%n_event=event_id
           if (event_id.gt.sim_obj(sid)%n_alloc_event) then
              sim_obj(sid)%n_alloc_event=&
                   sim_obj(sid)%n_alloc_event+&
                   alloc_event_n_in_sim_obj
              sim_obj(sid)%event_table_list=>&
                   reallocInt(sim_obj(sid)%event_table_list,&
                   sim_obj(sid)%n_alloc_event)
           end if
           sim_obj(sid)%event_table_list(event_id)=eid

           ! set the object 2
           sim_obj(cur)%n_event=sim_obj(cur)%n_event+1
           if (sim_obj(cur)%n_event.gt.&
                sim_obj(cur)%n_alloc_event) then
              sim_obj(cur)%n_alloc_event=&
                   sim_obj(cur)%n_alloc_event+&
                   alloc_event_n_in_sim_obj
              sim_obj(cur)%event_table_list=>&
                   reallocInt(sim_obj(cur)%event_table_list,&
                   sim_obj(cur)%n_alloc_event)
           end if
           sim_obj(cur)%event_table_list(sim_obj(cur)%n_event)=eid

           ! Set the event table
           event_table(eid)%obj1=sid
           event_table(eid)%obj_seq1=event_id
           event_table(eid)%obj2=cur
           event_table(eid)%obj_seq2=sim_obj(cur)%n_event
           event_table(eid)%event_type=rctn_event
           event_table(eid)%sqnc=k
           event_table(eid)%react_index=rid
           event_table(eid)%removed=.false.
           event_table(eid)%is_spontaneous=.false.

           ! Set the event profile
           if (react_type(rid)%react_const_rate(k).eq.0) then
              ! if the event is spontaneous,
              ! add to the spont_event_table
              event_table(eid)%is_spontaneous=.true.
              n_spont_event_table=n_spont_event_table+1
              spont_event_table(n_spont_event_table)=eid
              event_table(eid)%spont_table_index=&
                   n_spont_event_table
           else
              ! if the event is probabilistic,
              ! add to the event_profile
              if (profile_map_rctn(rid,k).eq.0) then
                 ! Create a new profile
                 n_event_profile=n_event_profile+1
                 call alloc_profile_table(n_event_profile,pid)
                 profile_map_rctn(rid,k)=pid
                 event_profile(pid)%event_type=rctn_event
                 event_profile(pid)%rate_const=&
                      react_type(rid)%react_const_rate(k)
                 event_profile(pid)%index=rid
                 event_profile(pid)%sqnc=k
                 event_profile(pid)%n_event=1
                 if (event_profile(pid)%n_alloc_event.eq.0) then
                    event_profile(pid)%n_alloc_event=&
                         alloc_event_n_in_profile
                    event_profile(pid)%event_list=>&
                         reallocInt(event_profile(pid)%event_list,&
                         event_profile(pid)%n_alloc_event)
                 end if
                 event_profile(pid)%event_list(1)=eid
                 event_table(eid)%profile_table_index=pid
                 event_table(eid)%profile_table_seq=1
              else
                 ! Add to the exisiting profile
                 pid=profile_map_rctn(rid,k)
                 n=event_profile(pid)%n_event
                 n=n+1
                 event_profile(pid)%n_event=n
                 if (n.gt.event_profile(pid)%n_alloc_event) then
                    event_profile(pid)%n_alloc_event=&
                         event_profile(pid)%n_alloc_event+&
                         alloc_event_n_in_profile
                    event_profile(pid)%event_list=>&
                         reallocInt(event_profile(pid)%event_list,&
                         event_profile(pid)%n_alloc_event)
                 end if
                 event_profile(pid)%event_list(n)=eid
                 event_table(eid)%profile_table_index=pid
                 event_table(eid)%profile_table_seq=n
              end if
           end if
        end do
100     cur=sim_obj(cur)%next
     end do
  end do

  ! event type 5 : sink with geometric object
  do i=1,n_geom_obj

     ! Check if the reaction is defined in the data sheet
     gid=sink_map(i,lid)
     if (gid.eq.0) cycle

     ! Calculate the distance between the geom_obj and the object
     select case(geom_obj(i)%geom_type)
     case(1)
        ! point
        hx=px-geom_obj(i)%x1
        hy=py-geom_obj(i)%y1
        hz=pz-geom_obj(i)%z1
        dist2=hx*hx+hy*hy+hz*hz
     case(2)
        ! straight line
        npos=((geom_obj(i)%x1-px)*(geom_obj(i)%x1-geom_obj(i)%x2)&
             +(geom_obj(i)%y1-py)*(geom_obj(i)%y1-geom_obj(i)%y2)&
             +(geom_obj(i)%z1-pz)*(geom_obj(i)%z1-geom_obj(i)%z2))/&
             ((geom_obj(i)%x1-geom_obj(i)%x2)**2&
             +(geom_obj(i)%y1-geom_obj(i)%y2)**2&
             +(geom_obj(i)%z1-geom_obj(i)%z2)**2)
        if (npos.lt.0) then
           hx=geom_obj(i)%x1
           hy=geom_obj(i)%y1
           hz=geom_obj(i)%z1
        else if (npos.gt.1) then
           hx=geom_obj(i)%x2
           hy=geom_obj(i)%y2
           hz=geom_obj(i)%z2
        else
           hx=npos*(geom_obj(i)%x2-geom_obj(i)%x1)+geom_obj(i)%x1
           hy=npos*(geom_obj(i)%y2-geom_obj(i)%y1)+geom_obj(i)%y1
           hz=npos*(geom_obj(i)%z2-geom_obj(i)%z1)+geom_obj(i)%z1
        end if
        dist2=(px-hx)**2+(py-hy)**2+(pz-hz)**2
     case(3)
        dist2=(px-geom_obj(i)%x1)**2
     case(4)
        dist2=(px-geom_obj(i)%x1)**2
     case(5)
        dist2=(py-geom_obj(i)%y1)**2
     case(6)
        dist2=(py-geom_obj(i)%y1)**2
     case(7)
        dist2=(pz-geom_obj(i)%z1)**2
     case(8)
        dist2=(pz-geom_obj(i)%z1)**2
     end select

     ! loop over the each type of sink reaction
     do j=1,sink_type(gid)%n_data

        ! Compare the distance and capture radius
        if (dist2.gt.sink_type(gid)%react_rd2(j)) cycle

        ! Determine the index of the event table
        n_event_table=n_event_table+1
        call alloc_event_table(n_event_table,eid)

        ! set the object
        event_id=event_id+1
        sim_obj(sid)%n_event=event_id
        if (event_id.gt.sim_obj(sid)%n_alloc_event) then
           sim_obj(sid)%n_alloc_event=&
                sim_obj(sid)%n_alloc_event+&
                alloc_event_n_in_sim_obj
           sim_obj(sid)%event_table_list=>&
                reallocInt(sim_obj(sid)%event_table_list,&
                sim_obj(sid)%n_alloc_event)
        end if
        sim_obj(sid)%event_table_list(event_id)=eid

        ! set the event table
        event_table(eid)%obj1=sid
        event_table(eid)%obj_seq1=event_id
        event_table(eid)%obj2=i
        event_table(eid)%obj_seq2=j
        event_table(eid)%event_type=sink_event
        event_table(eid)%react_index=gid
        event_table(eid)%removed=.false.
        event_table(eid)%is_spontaneous=.false.

        ! set the event profile
        if (sink_type(gid)%react_const_rate(j).eq.0) then
           ! if the event is spontaneous,
           ! add to the spont_event_table
           event_table(eid)%is_spontaneous=.true.
           n_spont_event_table=n_spont_event_table+1
           spont_event_table(n_spont_event_table)=eid
           event_table(eid)%spont_table_index=&
                n_spont_event_table
        else
           ! if the event is probabilistic
           ! add to the event_profile
           if (profile_map_sink(gid,j).eq.0) then
              ! Create a new profile
              n_event_profile=n_event_profile+1
              call alloc_profile_table(n_event_profile,pid)
              profile_map_sink(gid,j)=pid
              event_profile(pid)%event_type=sink_event
              event_profile(pid)%rate_const=&
                   sink_type(gid)%react_const_rate(j)
              event_profile(pid)%index=gid
              event_profile(pid)%sqnc=j
              event_profile(pid)%n_event=1
              if (event_profile(pid)%n_alloc_event.eq.0) then
                 event_profile(pid)%n_alloc_event=&
                      alloc_event_n_in_profile
                 event_profile(pid)%event_list=>&
                      reallocInt(event_profile(pid)%event_list,&
                      event_profile(pid)%n_alloc_event)
              end if
              event_profile(pid)%event_list(1)=eid
              event_table(eid)%profile_table_index=pid
              event_table(eid)%profile_table_seq=1
           else
              ! Add to the exisiting profile
              pid=profile_map_sink(gid,j)
              n=event_profile(pid)%n_event
              n=n+1
              event_profile(pid)%n_event=n
              if (n.gt.event_profile(pid)%n_alloc_event) then
                 event_profile(pid)%n_alloc_event=&
                      event_profile(pid)%n_alloc_event+&
                      alloc_event_n_in_profile
                 event_profile(pid)%event_list=>&
                      reallocInt(event_profile(pid)%event_list,&
                      event_profile(pid)%n_alloc_event)
              end if
              event_profile(pid)%event_list(n)=eid
              event_table(eid)%profile_table_index=pid
              event_table(eid)%profile_table_seq=n
           end if
        end if
     end do
  end do

  ! Add to the corresponding virtual box
  if (vbox(bid)%n_obj.eq.0) then
     vbox(bid)%first=sid
     sim_obj(sid)%next=sid
     sim_obj(sid)%prev=sid
  else
     sim_obj(sid)%next=vbox(bid)%first
     sim_obj(vbox(bid)%first)%prev=sid
     vbox(bid)%first=sid
     sim_obj(sid)%prev=sid
  end if
  vbox(bid)%n_obj=vbox(bid)%n_obj+1

  ! Count the object
  nx=ceiling(px/ch_width(1))
  ny=ceiling(py/ch_width(2))
  nz=ceiling(pz/ch_width(3))

  sim_obj(sid)%box(1)=nx
  sim_obj(sid)%box(2)=ny
  sim_obj(sid)%box(3)=nz

  do i=1,obj_type(lid)%n_group_belong
     n=obj_type(lid)%group_list(i)
     if (count_obj(n)%group_type.lt.3) then
        count_obj(n)%nx(nx)=count_obj(n)%nx(nx)+1
        count_obj(n)%ny(ny)=count_obj(n)%ny(ny)+1
        count_obj(n)%nz(nz)=count_obj(n)%nz(nz)+1
        total_count_obj(n)=total_count_obj(n)+1
     else
        nid=obj_type(lid)%config(count_obj(n)%point_type)
        count_obj(n)%nx(nx)=count_obj(n)%nx(nx)+nid
        count_obj(n)%ny(ny)=count_obj(n)%ny(ny)+nid
        count_obj(n)%nz(nz)=count_obj(n)%nx(nz)+nid
        total_count_obj(n)=total_count_obj(n)+nid
     end if
  end do

  return

end subroutine create_object
